Muscle RNA-seq analysis
Dmitry Ryazantsev
2026-01-02
1 RNA-seq of muscles
We use files from SYSTEMS BIOLOGY NII, IMBP, 2nd RNA-seq, ok quality.
(First RNA-seq failed because of EuroMint polymerase usage.) Fastq files
there:
/data7a/bio/human_genomics/muscle/imbp_muscles_generozov/
During sequencing, the samples were divided into four lanes to avoid cell fallout. This would result in bias, with rare transcripts underrepresented and common ones overrepresented. We don’t analyze this bias. Coverage is OK, and they also include exome and methylation information. Reads was not trimmed. Fastq read from 4 lanes was concatenated in 1 total file and quantified by Salmon.
metadata <- read.csv("/data7a/bio/human_genomics/muscle/imbp_muscles_generozov/meta_imbp_muscles.csv")
tx2gene <- read.csv("/data7a/bio/human_genomics/muscle/analysis/timur_files/itog/tx2gene.csv")
met <- read.csv("/data7a/bio/human_genomics/muscle/analysis/timur_files/itog/samples.csv", stringsAsFactors = F)
met$Condition <- as.factor(met$Condition)
met$Person <- as.character(met$Person)
rownames(met) <- met$Sample
file_n <- setNames(met$Filepath, met$Sample)Expression matrix
txi <- tximport(
file_n,
type = "salmon",
tx2gene = tx2gene,
ignoreTxVersion = TRUE,
countsFromAbundance = "no")
saveRDS(txi, '/data7a/bio/human_genomics/muscle/analysis/RDS/txi.RDS')Lets check values’ distribution. Попробуем разделить денсити по образцам. We observe matrix in base state, not in log state.
2 PCA plots
Lets plot pca to check bias.
log_pca <- log10(tpm_matrix + 1)
pca <- pca(log_pca, metadata = met, removeVar = 0.1)screeplot(pca, components = getComponents(pca, seq_len(6)))
col_by <- "Condition"
biplot(
pca, x = "PC1", y = "PC2", pointSize = 3, colby = col_by,
legendPosition="right", showLoadings = F)
Findings from PCA:
- There is an artifact sample 9_B1, it needs to be counted with and without it.
- You can try to identify the tissue of 9_B1.
- Points B1 and B2 are divided into groups, but half of the B2 points are mixed into the B1 group. This indicates either strong individual differences (e.g., different muscle growth rates—this requires calibrating the measurement instruments) or a technical error (missing the target and taking the wrong tissue, poor rRNA depletion, insufficient sample size).
3 Differential Expression
Lets count with and without 9B outliers samples.
dds <- DESeqDataSetFromTximport(
txi,
colData = met,
design = ~ Person + Condition)
keep <- rowSums(counts(dds) >= 10) >= ncol(dds)/2
dds <- dds[keep, ]
dds <- DESeq(dds)
res_all <- results(dds, contrast = c("Condition", "Trained", "Untrained"))
# Фильтрация значимых генов
sig_genes_all <- as.data.frame(res_all)[!is.na(res_all$padj) & res_all$padj < 0.05 & abs(res_all$log2FoldChange) > 1, ]
# add symbols
sig_genes_all$symbol <- mapIds(org.Hs.eg.db, keys = rownames(sig_genes_all),
column = "SYMBOL", keytype = "ENSEMBL", multiVals = "first")
sig_genes_all <- sig_genes_all[order(sig_genes_all$padj), ]Try without weird sample.
met_filt <- met[met$Sample != "9_B1", ]
samples_to_keep <- met_filt$Sample
txi_filt <- txi
txi_filt$counts <- txi$counts[, samples_to_keep]
txi_filt$abundance <- txi$abundance[, samples_to_keep]
txi_filt$length <- txi$length[, samples_to_keep]
rm(samples_to_keep)
dds2 <- DESeqDataSetFromTximport(
txi_filt,
colData = met_filt,
design = ~ Person + Condition)
keep <- rowSums(counts(dds2) >= 10) >= ncol(dds2)/2
dds2 <- dds2[keep, ]
dds2 <- DESeq(dds2, parallel = T)
res_nooutl <- results(dds2, contrast = c("Condition", "Trained", "Untrained"))
# Фильтрация значимых генов
sig_genes_nooutl <- as.data.frame(res_nooutl)[!is.na(res_nooutl$padj) & res_nooutl$padj < 0.05 & abs(res_nooutl$log2FoldChange) > 1, ]
# add symbols
sig_genes_nooutl$symbol <- mapIds(org.Hs.eg.db, keys = rownames(sig_genes_nooutl),
column = "SYMBOL", keytype = "ENSEMBL", multiVals = "first")
sig_genes_nooutl <- sig_genes_nooutl[order(sig_genes_nooutl$padj), ]
write_csv(sig_genes_nooutl, '/data7a/bio/human_genomics/muscle/analysis/download/DEGs_nooutl.csv')Check intersection between outliers filtered and unfiltered.
v <- list(Filtered = rownames(sig_genes_nooutl), `All samples` = rownames(sig_genes_all))
ggVennDiagram(v) + scale_fill_gradient(low = "#F4FAFE", high = "#4981BF") + coord_flip()
3.1 Compare two ways
print("Without outlier")## [1] "Without outlier"summary(res_nooutl)##
## out of 14528 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 979, 6.7%
## LFC < 0 (down) : 592, 4.1%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsprint("All samples")## [1] "All samples"summary(res_all)##
## out of 14582 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 868, 6%
## LFC < 0 (down) : 471, 3.2%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsThe DEGs number increased, not decreased. If the outlier carried a biological signal, there would have been fewer of them. Here, it’s the opposite—the outlier smeared the variance, weakened statistical power, and obscured the effects between groups.
MA-plot — is a diagram in which each point is a gene.
A (average expression) — average gene expression level (log scale)
M (minus / fold change) — logarithm of the expression ratio between two conditions (log2FoldChange)
A = 0.5 * log2(baseMean) M = log2FoldChange(sample2 / sample1)
plotMA(res_nooutl, ylim = c(-5, 5), main = "MA-plot Filtered")
plotMA(res_all, ylim = c(-5, 5), main = "MA-plot Unfiltered")
We need to look at the cloud’s displacement to see if the normalization is correct.
Conclusion from the MA plot: the Filtered values are more closely clustered, meaning the dispersion has decreased, meaning the outlier was indeed noisy.
Dispersion estimates for filtered and unfilt.
plotDispEsts(dds, main = "Dispersion estimates unf")
plotDispEsts(dds2, main = "Dispersion estimates flt")